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ABSTRACT 

Weak-lensing searches for galaxy clusters are plagued by low completeness and 
purity, severely limiting their usefulness for constraining cosmological parameters 
with the cluster mass function. A significant fraction of 'false positives' are due to 
projection of large-scale structure and as such carry information about the matter 
distribution. We demonstrate that by constructing a "peak function" , in analogy 
to the cluster mass function, cosmological parameters can be constrained. To this 
end we carried out a large number of cosmological A^-body simulations in the 
^m-cs plane to study the variation of this peak function. We demonstrate that 
the peak statistics is able to provide constraints competitive with those obtained 
from cosmic-shear tomography from the same data set. By taking the full cross- 
covariance between the peak statistics and cosmic shear into account, we show 
that the combination of both methods leads to tighter constraints than either 
method alone can provide. 
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1 INTRODUCTION 

The number density of clusters of galaxies is a sen- 
sitive probe for the total matter density of the Uni- 
verse f2 m , the normalisation of the power spectrum as, 
and the evoluti on of the equation of sta t e of the Dark 
Energy w (e.g.. IWang fc SteinhardtJll998t H aiman et al.l 
l200ll ; IWeller et alj|2002l ). For some time it was thought 
that weak gravitational lensing, which by its nature is 
sensitive to dark and baryonic matter alike and inde- 
pendent of the dynamical or evolutionary state of the 
cluster, could be used to construct clean, purely mass 
selected cluster samples. However, ray-tracing simula- 
tions through cosmological Af-body simulation made it 
clear that weak-lensing selected clusters are not at all 
mass selected but selected by the she ar of the pro- 
jected mass along the line of sight (e.g., Hamana et al.l 
12004 iHennawi fc Spergell 120051 ; IDietrich et al.ll2007l ). As 
a result, blind searches for galaxy clusters using weak 
lensing have both low purity and com pleteness (e.g., 
ISchirmer et~aL1l2007l ; IDietrich et "aL1l2007l 1. 

Gravitational lensing is, due to the large intrin- 
sic ellipticity scatter of background galaxies, an inher- 
ently noisy technique. This shape noise is the domi- 
nant noise source at the low signal-to-noise ratio (SNR) 
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end of the weak-lensing selection function, while pro- 
jections of large-scale structure (LSS) along the line-of- 
sight (LOS) dominate the noise budget of highly signifi- 
cantly detected peaks (|Dietrich et al.l2007t ). Both sources 
of noise affect purity and completeness. Galaxy clusters 
aligned with underdense regions are not visible as sig- 
nificant overdensities, while the projection of uncorre- 
cted overdensities can mimic the shear signal of galaxy 
clusters. While these effec ts can be taken into account 
dMarian fc Bernstein 2006), they degrade the constraints 
on cosmological parameters one can obtain using weak- 
lensing selected galaxy clusters. 

Of course such projected peaks are noise or false 
positives only in the sense of galaxy cluster searches. 
They are caused by real structures along the line-of-sight 
and as such carry information about the matter power 
spectrum. Wh ereas analytical models exist for the halo 
mass function (|Press fc Schechterlll974l ; ISheth fc Tormenl 
2002), no such model exists for the number density of 
peaks in weak lensing surveys. Probably no such predic- 
tion can be made analytically because the abundance of 
peaks depends on projections of uncollapsed yet highly 
non-linear structures like filaments of the cosmic web. 
As an additional complication the observed number of 
peaks depends on observational parameters like limiting 
magnitude, redshift distribution, and intrinsic ellipticity 
dispersion. 

In the absence of an analytic framework, ray-tracing 
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through A-body simulations can be used to numerically 
compute the "peak function" (in analogy to the mass 
function) for a survey and study its variation with cos- 
mological parameters. Here we present a large set of 
such simulations aimed at demonstrating the usefulness 
of the shear-peak statistics for constraining cosmolog- 
ical parameters. We consider this work to be a pilot 
study and limit ourselves to the variation of the peak 
function with f2 m and as and its ability to break the 
degeneracy between these two parameters encountered 
in the 2-poin t cosm ic-shear correlation-function. Unlike 
iMarian et all (2009) who showed that the projected mass 
function, which is difficult to measure, scales with cosmol- 
ogy essentially in the same way as the halo mass function, 
we study the cosmological dependence of the directly ob- 
servable aperture mass statistics. 



2 METHODS 

2.1 A-body simulations 

We carried out A-body simulations for 158 different flat 
ACDM cosmologies with varying S7 m , S7a, and erg. Fig- 
ure [1] shows the distribution of these simulations in the 
Om-as plane. All simulations had 256 3 dark matter parti- 
cles in a box with 200 hy Mpc side length. These choices 
reflect a compromise we had to make between computing 
a large number of simulations to sample our parameter 
space on the one hand and to have a fair representation of 
very massive galaxy clusters dominating the cosmological 
sensitivity of the halo mass function on the other hand. 
These simulation parameters were chosen such that we 
can expect the presence of 10 15 hjg Mq mass halos at 
redshift z — in the simulation box in our choice of 
fiducial cosmology 7T = (fi mo — 0.27, Qa = 0.73, fib = 
0.04, cr 8o = 0.78, n s = 1.0, T = 0.21, h 70 = 1). 

We computed 35 A-body simulations for this fiducial 
cosmology to estimate the covariance of our observables. 
The total number of A-body simulations is thus 192. Par- 
ticle masses depend on the cosmology and range from 



9.3 x 10 a M Q for Q. 



0.07 to m D 



,2 x 10 10 M, 







for Q m — 0.62. The particle mass at our fiducial cosmol- 
ogy is m p = 3.6 x 10 10 Mq. 

The A-body simulations were carried out with the 
publi cly available TreePM code GADGET-2 (|Springell 
20051). The initial con ditions were generated using the 
Eisenstein fc Hul (|l998T ) transfer function. We started the 
simulations at z — 50 and saved snapshots in Az inter- 
vals corresponding to integer multiples of the box size, 
such that we have a suitable snapshot for each lens plane 
of the ray-tracing algorithm. The Plummer-equivalent 
force softening length was set to 25 h^g kpc comoving. 
We checked the accuracy of our A-body simulations by 
comparing their matte r powe r spectra with the fitting 
formula of ISmith et al.l (|2003l ) . Additionally, we also de- 
tected halos using a friend-of-friend halo finder and com - 
pared their mass function to that of I Jenkins et al.l l|200ll ') . 
All tests were done for a number of different cosmologies 
and redshifts to ensure that the simulations match our 
expectations over the parameter and redshift range un- 
der investigation here. 
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Figure 1. Location of the 158 different cosmologies in the 
fim-o"8 plane for which iV-body simulations were computed. 
The red diamond marks the fiducial cosmology at (Qm, <rg) = 
(0.27,0.78). 



2.2 Ray-tracing 

We used the mult iple lens -plane algorithm (e.g . 

198fj ; ISchneider et~aLl 1 19921 ; 

[20 



Blandford & Naravan 



Seitz et all 1 1994 I Jain et al.l l2000l ; IHilbert et all l2009h 
to simulate the propagation of light rays through the 
matter distribution provided by the A-body simulations: 
for a given A-body simulation, we constructed the 
matter distribution along the line of sight by tiling 
snapshots of increasing redshift. The matter distribution 
of each snapshot was projected onto a lens plane located 
at the snapshot redshift. 

Note that the boxes are just small enough for the 
cosmic evolution during the light travel time through a 
box to be negligible. This ensures that the matter distri- 
bution does not change significantly in the volume that 
is represented by a particular lens plane, and that the 
scale factor and the comoving angular diameter distances 
to the structure projected onto this plane are essentially 
the same. If the latter were not the case, this would lead 
to an erroneous conversion of physical scales on the lens 
plane to angular scales on the sky. Our allowed us to 
project a complete snapshot onto one lens plane instead 
of creating several smaller redshift slices as was done in 
IHilbert et ail (2009). which reduces the complexity of the 
ray-tracing considerably. 

Since the snapshots basically contain the same mat- 
ter distribution at slightly different stages of evolution, 
measures have to be taken to avoid the repetition of struc- 
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tures along the line of sight. Making use of the periodic 
boundary conditions of the simulation volume, we ap- 
plied random rotations, translations, and parity flips to 
the matter distribution of each snapshot prior to the pro- 
jection. For the ray-tracing, we assumed that light rays 
are only deflected at the lens planes and propagate freely 
in between. We compute the Fourier transform of the de- 
flection potential on each plane from the projected mass 
density by solving the Poisson equation in Fourier space 
using FFT, again exploiting the periodic boundary condi- 
tions. From this, the Fourier transforms of the deflection 
angles and their derivatives can be obtained using simple 
multiplications. Finally, these quantities are transformed 
to real space using an inverse F FT. More details on the 
formalism can be found, e.g., in IJain et all l|2000l ). With 
this, a set of light rays (forming a grid in the image plane) 
can be propagated from the observer throu gh the array of 
lens planes using a recursion formula (see iHilbert et al.l 
2009). Similarly, the Jacobian matrix of the lens map- 
ping from the observer to each of the lens planes can be 
obtained. 

We then sampled the image plane uniformly with 
galaxies, the redshift of which was drawn from a distri- 
bution of the form 



0) « 



exp 



(1) 



The Jacobian matrices were interpolated from the grid 
onto the galaxies (in the plane of the sky as well as in 
redshift) and the reduced shear was computed. We sim- 
ulated a CFHTLS-Wide like survey for which we created 
five 6 x 6 sq. deg. patches from every TV-body simulation. 
The parameters of the redshift distribution ([!} were set 
to a = 0.836, /? = 3.425, and zp = 1.171, as de termined 
for the CFHTLS-Wide ijBeniamin et al.l 12007^ . We set 
the galaxy number-density to n g = 25 arcmin -2 and the 
intrinsic ellipticity dispersion to o> = 0.38. Because very 
few galaxies are present at high redshifts, the redshift 
distribution was cut off at z — 3.0 to save computing 
time. 



2.3 Peak detection 

2.3.1 Aperture mass in 2-d 

The tidal gravitational field of matter along the line- 
of-sight causes the shear field 7(0) to be tangentially 
aligned around projected mass-density peaks. We can 
use this tangential alignment directly to detect weak- 
lensing peaks, instead of searching for convergence peaks 
on maps of reconstructed surface mass-density as it has 
been done sometimes in weak-lensing cluster sea rches 
(e.g.. lGavazzi fc Soucau1l2007l; iMivazaki et afll2007l ). We 
define the aperture mass (|Schneiderlll996l ) at position Go 
to be the weighted integral 



>(*<>) = J 



d 2 0QO?ht(0;0o) 



(2) 



supQ 



over the shear component tangential to the line Go — 
0, 7t(0;0 o ). Here Q{&) = Q(\0\) is a radially sym- 
metric, finite and continuous weighting function with 



limtf-.oo Q(i?) = 0. For later convenience we also require 
Q to be normalised to unit area. If Q(#) follows the ex- 
pected shear profile of a mass peak, the aperture mass 
becomes a matched filter technique for detecting such 
mass peaks. On data the shear field is sampled by galax- 
ies with ellipticities e%. Then M ap can be estimated by 
the sum over N s galaxies in the aperture, 



Map = -V<3(l),)£,t 

TV „ ^ ' 



(3) 



where e« is the tangential ellipticity component of the 
i-th galaxy, defined in analogy to jt above. 

The SNR of the aperture mass can be computed di- 
rectly from the data, making use of the fact (Map) = 0. 
Then the RMS dispersion is o"M ap = \J (A/ ap ) 2 , which can 
be estimated by 



V2r, 



1/2 



where we have made use of the fact that 

2 



(4) 



(•») 



with o> being the intrinsic ellipticity dispersion. The es- 
timator for the SNR of the aperture mass is then finally 



5(0o ) 



(6) 



2.3.2 Tomographic aperture mass 

The aperture mass statistics locates convergence peaks 
only in projection on the sky. Using redshift information 
on the background galaxies, e.g., from photometric red- 
shifts, one can generalise the 2-dimensional aperture mass 
to a tomographic measure that is able to deproject struc- 
tures along the line-of-sight and locate peaks in redshift 
space (IHennawi fc Spergel|[2005h . The likelihood that a 
peak at a position Go is at a redshift z& is given by 

ln£(0 O) * d ) = ^ E „. Z2(zi;zd) , (7) 

where Z(zi\Zd) is the redshift weight for a background 
galaxy in the ith redshift bin, 



Z(z; 2 d ) = Dd ^ ds H(z - z d ) 



(8) 



with the Heaviside step function H. A peak is then lo- 
cated at the 3-d position (0o, Zd) that maximises the like- 
lihood £.. For the purpose of this work 10 equally spaced 
redshift steps Zd ~ 0.1 . . . 1.0 were used. The background 
galaxies were put into redshift bins with width Az = 0.01 
assuming perfect knowledge of their redshifts. 

As in (Dietrich et all (120071) we used the weight func- 
tion proposed bv lSchirmer et al.l (|2007l ) 



Jnfw(i;ic) oc 



ta,nh(x/x c 



1 + e 6 " 150 * + e - 47 + 5to x/x c 



(9) 

where x — $/i9 max and x c is a free parameter, which was 
fixed to the value of x c — 0.15 determined to be ideal for 
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the d etection of galaxy clusters by iHetterscheidt et alj 
(2005). Qnfw follows the shear profile of an NFW halo 
with exponential cut-offs as x — > or x — > oo. 

The absolute scale i? max determines the halo radius 
or mass to which the filter function is tuned. The filter 
scale chosen for our simulations is 5'. 6 on the sky, corre- 
sponding to a radius of Mpc at a redshift of z — 0.3. 
At this redshift the lensing efficiency of our survey is 
maximal and the chosen radius is adjusted to cluster 
sizes easily detectable with weak lensing while smoothing 
over smaller halos. This smoothing also ensures that shot 
noise from unresolved structures in the JV-body simula- 
tions does not play a role. 

Peaks were detected by connected-component la- 
belling of pixels above a detection threshold. We used 
the 8 connectivity in 2-d and the 26 connectivity in 3- 
d, i.e., we consider all pixels that are connected via the 
sides, edges, or corners of a square or a cube as one struc- 
ture. Additionally, for tomographic peaks the condition 
was imposed that peaks must be detected in at least 
three adjacent redshift bins. This additional requirement 
is used to filter out detections at very high or low redshifts 
whose true redshift is outside the tomography cube. Such 
peaks would pile up in the lowest and highest redshift bin 
and lead to high additional noise in them. At the same 
time this filter criterion suppresses the inclusion of peaks 
caused by increasing shot noise at high redshifts caused 
by the sharply decreasing number density of background 
galaxies. Such peaks occur typically in only one or two 
redshift slices. 

2.4 Analysis 

For every cosmological model the peak function gives a 
p-dimensional data vector £ of observables. We will ex- 
plore several choices of observables below. The choice of 
cosmological parameters is denoted by 7r and the model 
prediction is m(n). The posterior probability distribu- 
tion is 



P(n\0 = ^gW) 



(10) 



where p(7r) is the prior probability distribution, p(CI 7r ) 
is the likelihood, and p(C) is the evidence. We used a flat 
prior with cutoffs, i.e., p(n) = 1 if fi m £ [0.1 : 0.5] and 
erg G [0.4 : 1.1] and p(n) = otherwise. The evidence 
in our case simply is a normalisation of the posterior ob- 
tained by integrating the likelihood over the support of 
the prior. 

Assuming that £ has a Gaussian distribution, the 
likelihood is 

p«k) = 1 

(2iv) d / Vdct E(tt) 

x exp |-i [C - m(TT)] t T' 1 (-n-) [£ - m(n) 

(11) 

where E(tt) is the covariance matrix of the d-dimensional 
vector £. Since our parameter space is covered only by 
discrete points we will compute £ by fitting smooth func- 
tions to our data vectors or by interpolating across our 
parameter space. Details will be given in Sect. [3] 



We obtained estimates C of the covariance matri- 
ces for the shear tomography and the peak statistics, as 
well as their cross-covariance, from the field-to-field vari- 
ation in the 175 ray-tracing simulations for the fiducial 
cosmology. As indicated in Eq. pip , the covariance in 
principle depends on cosmology. Since we do not have 
a sufficient number of simulations for other cosmological 
parameters, we set E(7r) = E(7To). Although this is an 
approximation commonly made, neglecting the cosmol- 
ogy dependence of E can have a non-negligible impact on 
the shap e of the posterior l ikelihood, as has been investi- 
gated in lEifler et al.l (|2008h for the case of cosmic shear. 
Furthermore, we note that the assumpti on of a Gaus- 
sian l ikelihood is not necessarily justified ()Hartlap et al.l 
2009; Schneider fc Hartladliooa) . These studies suggest 
that both approximations lead to an over-estimation of 
the errors on the cosmological parameters, so that our 
confidence regions are most likely very conservative. 

A further issue to keep in mind when estimating a 
covariance matrix from the data or from simulations is 
that, although the covariance C estimated from the data 
is an unbiased estimator of the true covariance E, the 
inverse C _1 is not an unbiased estimator of E _1 . For n 
independent simula tions an unbiased es timator of the in- 
verse covariance is (jHartlap et alJl2007l ) 



c- 



d- 



n-1 



2 c- 



(12) 



which is what we used when computing Eq. pip . How- 
ever, the estimated covariance matrix becomes singular 
if d > n— 1, which means that the limited number of ray- 
tracing simulations available to us constrains the number 
of bins that can be used for the analysis. 



3 RESULTS 

In analogy to the mass function N(M, z\n), the peak 
function measures the abundance of peaks as a function 
of convergence and redshift N(k, z\ir), where for single 
structures along the line-of-sight the variation of the con- 
version from k to M with redshift is given by the ker- 
nel (|8]). Because we detected peaks not in convergence 
maps but in aperture-mass maps or cubes, the SNR of 
peaks was used as a proxy for mass. 



3.1 Constraints from aperture mass maps 

As a first step we show that constraints on 7r = (f2 m , (Ts) 
can be obtained from the peak statistics in the absence 
of redshift information. A sim ilar study was recently pre- 
sented by IWang et al.l ([2009), who demonstrated that 
parameters of the Dark Energy equation of state can 
be constrained from high convergence regions. We de- 
tected peaks in aperture-mass SNR maps as described in 
Sect. 12.41 with a detection threshold of 3.25<r. 

Binning the peaks by SNR is not an ideal way to 
use information about their projected mass since either 
high SNR bins in cosmologies with low clustering remain 
empty, or very massive peaks are beyond the SNR of the 
maximum bin. Instead we used the cumulative SNR dis- 
tribution of peaks. The function S(n m ,as) : R 2 -> l" bi ° 
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SNR 

Figure 2. Construction of the function S. The solid black line 
is the cumulative SNR distribution of peaks detected in one 
of our 35 realizations of the fiducial cosmology. The horizontal 
dashed lines are the logarithmically spaced percentiles from 
/min = 0-5 to / max = 0.98 at which the cumulative SNR dis- 
tribution is sampled. The corresponding SNR values denoted 
by the vertical dashed lines are the values in our data vector. 



gives the SNR at which the cumulative distribution ex- 
ceeds the /th percentile for 7ibi n values of / ranging from 
/min to /max- Figure [2] illustrates how S is constructed. 

We measured S(fl m ,os) for n hin — 5 logarithmically 
spaced values from / m i n = 0.50 to f max = 0.98. At the 
fiducial cosmology these percentiles corresponds to SNR 
values of 3.5cr and 5.7a, respectively. Typically several 
hundred peaks per 36 sq. deg. field were detected so that 
the 98%ile could be reliably measured. 

We used bilinear smoothing splines (|Dierckxlll993l ) 
to interpolate S(Q m ,ffs) on the grid covered by our N- 
body simulations. In this section splines are a sufficient 
description of the variation of S over our parameter space 
because we only seek to qualitatively demonstrate the 
ability of the peak statistics to constrain cosmological 
parameters and to illustrate some of its properties. We 
will use a more quantitative approach in the following 
sections. 

Figure [3] shows the confidence contours derived from 
this statistics in the S7 m -o"s plane. They have a shape 
similar to that see n in constraints d erived from clus- 
ter cosmo l ogy (e.g., iHenrv et alJ feuO^ and cosmic shear 
(|Fu et all [200a . e.g.,) for a CFHTLS like 180 sq. deg. 
surey. In order to achieve this, we scaled the covariance, 
which we computed for the individual 36 sq. deg. fields 
back to the full survey. The similarity of the constraints 
is of course no surprise since the peak statistics measures 
the same density fluctuations as clusters of galaxies and 
cosmic shear. 

Although the spline interpolation is mostly illustra- 
tive, we defined a figure of merit (FoM) , in analogy to 
the F oM of the Dark Energy Task Force (|Albrecht et al.1 
2006), as the inverse of the area inside the 95% confidence 
contour. We used this FoM to characterize how the peak 
statistics changes when parameters entering the function 
S are modified. Here in particular we examined the de- 
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Figure 3. Confidence contours of the aperture mass peak 
statistics. Shown are the 1-, 2-, and 3<r confidence contours 
of the S statistics. The white cross denotes the fiducial cos- 
mology. 



pendence of the cosmological constrains on the minimum 
significance of a detection. 

The detection threshold employed in the produc- 
tion of Fig. [3] is very low and a sizable fraction of the 
peaks detected in thi s way are simply due to shape noise 
l|Dietrich et al.l [20071 ) and do not carry cosmological in- 
formation. However, at such a low detection threshold 
most peaks not caused by noise fluctuations are also not 
due to a single massive halo but caused by the alignment 
of LSS along the LOS. We demonstrate that these low 
significance peaks indeed carry cosmological information 
by comparing the FoM of the statistics in Fig. [3] to the 
FoM resulting from the same function S with a detection 
threshold of 4.5a. While the constraints in Fig. [3] cor- 
respond to a FoM of 40, the higher detection threshold 
results in a FoM of only 20. We note that the 95% con- 
fidence interval is not fully contained in the support of 
our flat prior. For the low SNR detection, the 95% con- 
fidence interval is cut off by the prior only at the high 
f2m/low as end. The prior terminates the banana shaped 
confidence region at both ends for the high SNR detec- 
tion constraints. Consequently, the true figures of merit 
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for these statistics are smaller than presented here but 
more so for the higher detection threshold, making the 
importance of low SNR peaks even more evident. How- 
ever, because we used a very simple interpolation scheme, 
these numbers can only be rough guidelines and we will 
present a more detailed assessment of the S function in 
the next section. 

We emphasize that, although there is no reason to 
believe that the true values of O m or a$ are outside the 
support of our prior, the prior is used only to avoid ex- 
trapolating beyond the parameter range covered by the 
iV-body simulations. Since the aim of this study is to 
examine how well the peak statistics can constrain cos- 
mological parameters, we did not regard the prior as in- 
formation that should be included in the calculation of 
the FoM. 



3.2 Constraints from peak tomography 

In this section we make use of the redshift information 
in our shear catalogues. We employed the tomographic 
aperture mass outlined in Sect. 12.331 to locate high con- 
vergence regions not only in projection on the sky but 
also along the redshift axis. With 175 independent ray- 
tracing simulations we could not compute the covariance 
of the full peak function N(S, z) for meaningful num- 
ber of bins in SNR and redshift. Instead, we constructed 
two separate peak functions from the tomographic data 
cubes. 

The first function measures the abundance of peaks 
in every redshift bin as a function of cosmology only. We 
detected peaks as described in Sect. 12.3.21 with a min- 
imum detection threshold of a m i n = 3.2 and assigned 
them to a redshift bin based on the redshift z& that 
maximises the likelihood (|7]l. With 10 redshift bins, the 
vector-valued function M(fl m , a&) : R 2 — > R 10 counts 
the number of peaks in each redshift bin as a function of 
cosmological parameters. The second function S uses the 
SNR information of the detected peaks as defined in the 
previous section. We used the same values for / m in, /max, 
and Ubin. 

To interpolate M(Cl m ,ffs) and S , (r2 m ,o"8) between 
points covered by our iV-body simulations, we deter- 
mined fitting functions. These are given in Appendix |A1 
and are typically accurate to w 10%. Figure [4] shows con- 
fidence contours in the Q m -as plane obtained using these 
fitting functions individually and for the combination of 
both peak functions, £ p = (M , S) 1 . As one would expect 
from cosmic shear and the cluster mass function, there 
is a significant degeneracy between f2 m and <T8. The in- 
teresting result in Fig. [4] is that M and S 1 contain com- 
plementary information such that their degeneracies are 
partly broken. Although the kernel |(HJ is very broad and 
the deter mination of a peak's redsh ift consequently noisy 
(see also Irlennawi fc Spergell I2005T ) the information on 
the growth of structure with redshift can be statistically 
recovered with peak tomography. This then provides in- 
formation not contained in S, which is equivalent to a 
projected mass function. It is important to emphasise 
that we utilised the full cross-covariance between M and 
S when we computed the joint confidence contours. 




Figure 4. Confidence contours of the peak statistics. Green 
and blue contours are 1- and 2a contours for the M and S 
statistics, respectively. Limited by black contour lines are the 
joint 1-, 2-, and 3<r confidence intervals. 



We now revisited the issue of how the various param- 
eters of the S function change the information content of 
the peak statistics. In the previous section we showed that 
the detection threshold is an important parameter and 
that the number of projected low SNR peaks helps to con- 
strain parameters. The same is true for the tomographic 
peak statistics. The FoM of the S-function in Fig.|4]is 48. 
Increasing the detection threshold to 4cr and keeping all 
other parameters unchanged decreases the FoM to only 
7. Like in the case of the aperture-mass peak statistics, 
the confidence region is terminated by the prior and the 
true difference between the different detection thresholds 
is even larger than suggested by the FoM. 

We note that the FoM of the S function for peak 
catalogues generated from tomography cubes is not much 
higher than the rough estimate of the FoM in Sect. 13. ll of 
the projected peak statistics. The 95% confidence inter- 
val in the tomographic case is not fully contained in the 
parameter space explored by our simulations, whereas, 
with the spline approximation from the previous section, 
the low f2 m /high ag end of the confidence contours is 
within our parameter range. In any case, the substantial 



© 2009 RAS, MNRAS 000, QMT] 



Cosmology with the shear-peak statistics 7 



Table 1. Figures of Merit of different parameters of the S 
function. 



Table 2. Cosmological constraints using different statistics. 



/min 


/max 


"bin 




FoM 


0.50 


0.98 


5 


3.2 


48 


0.50 


0.98 


5 


4.0 


7 


0.08 


0.98 


10 


.3.2 


48 


0.08 


0.50 


5 


.3.2 


34 



gain of the tomographic peak statistics does not come 
from the deprojection of structures along the LOS but 
from localising peaks along the redshift axis, i.e., from 
the combination of the S and M statistics. 

We also tested several choices of / m i n and / max and 
the number of bins; Table [1] gives an overview of various 
settings. Extending / m i n to lower values adds almost no 
information, even if the number of bins is increased to 
preserve the information in the high SNR bins. For ex- 
ample, decreasing / m i n = 0.08 and setting nt,m = 10 does 
not improve the FoM. The information content of the S 1 
function is slightly higher with these settings, as is evi- 
denced by a 3% decrease of the area inside the 68% con- 
fidence interval. However, the 95% confidence contours 
within the support of our prior are not tightened. 

Despite of this, most of the information is contained 
in the low SNR regime, as can be seen from the last two 
rows of Table [TJ This information, however, can be re- 
covered with only one or two bins: Figure [2] shows that 
the low SNR end of the cumulative distribution function 
has an almost constant slope. This slope is completely 
determined by the SNR sampled at f m i n and the point 
(o"min,0), and explains why adding more points at the 
low end does not increase the FoM. Further information 
comes only from the shape of the curve in the interme- 
diate SNR range. At the high SNR end the cumulative 
distribution again has a constant slope and is fully char- 
acterized by the last two sample points. 



3.3 Combination with cosmic shear 

Up to now cosmological information has generally been 
extracted from lensing surveys via the c osmic-shear two - 
point correlation functions (2PCF) fe.g.,lSclmcidcr l2006l ). 

i±(9) = <et(0)et(«> + *)) ± (ex {0)e x (6 + ■&)) , (13) 

which ar e related to t he convergence power spectrum V K 
via (e.g., lKaiserlll992T ) 



r°° di 



dZ 
2t7 



Ji(W)V K (l) 



(14) 
(15) 



Here e x is the cross-component of the ellipticity and 
the J n (x) are the n-th Bessel functions of the first 
kind. Surveys using this method have led to in- 
creasingly tight c onstraints in the Qm - ag plane (e.g., 
Jarvis et al . 2006; Scmbol oni et alj|2006l; iHoekstra et al.l 
20061: iHetterscheidt et all 120071 ; iBeniamin et al.1 120071 : 



Fu et al.l l200Sf ) . However, the cosmic shear 2PCF de- 



Type 
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^8 
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Cosmic shear 
Peak statistics 
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776+ ' 107 
u -' '"-0.096 
774+0.095 
u - ' ' 4 -0.087 


71 
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pletely if they are purely Gaussian. Cosmic shear can ac- 
cess information about the non-Gaussianity of the matter 
distribution o nly through higher-o rder correlation func- 
tions (see e.g.. lTakada fc Jainll2003l . for constraints using 
the three-point correlation function). The peak statistics 
on the other hand is most sensitive to extreme overdensi- 
ties along the LOS, i.e., to those structures that contain 
most information about non-Gaussianity. It is thus rea- 
sonable to assume that both statistics are not completely 
degenerate and that combining the two does not sim- 
ply amount to using the same information twice. This 
expectation is supported by a number of studies look- 
ing at the constraints one can place on the Dark Energy 
equation of state by combining the cluster mass func- 
tion with other c osmological probes, including weak grav- 
itational lensin g (Fang & Haiman 2007; Takada ~fc Bridle] 
l2007l : ICunha et : alJl2009h . lTakada fc Bridle! (|2007l ) in par- 
ticular examined the full cross-covariance between the 
cosmic shear 2PCF and cluster counts of shear-selected 
halos and found that the combination of both methods 
leads to tighter constraints than either method alone can 
provide. 

Cosmic shear, like the peak statistics, can greatly 
benefit from the inclusion of redshift informati on by fol- 
lowing the evolut ion of structure with time l|Hul Il999l : 
iBacon et al.ll2005h . This is done by dividing the galaxy 
sample into redshift bins and computing their auto- and 
cross-correlation functions, 



scribes the underlying density fluctuations only com- 



CTV) = (e[ t) (0)e^(e + ^)±(e^(e)e^(e + ^)) . (16) 

We used this tomographic 2PCF to compare and combine 
the constraints obtained from cosmic shear to those from 
the peak statistics in the same survey. 

We split the galaxy catalogue into two redshift bins, 
separated at redshift z = 0.6 and estimated the to- 
mographic cosmic shear 2PCF C± 

in our simulation 

of the fiducial cosmology for 60 logarithmically spaced 
intervals from 30" to 6°. From this we constructed 
data vectors by averaging the values of 6 consecutive 
bins into one entry in the data vector, so that wc 
have a 60-dimensional cosmic shear data vector £ cs = 
($ 1 \£ 1 \$ 2 \e 2) ,i? 2 \e 2) T- Using these vec- 
tors we computed the covariance of our cosmic shear mea- 
surements at the fiducial cosmology. The choice of red- 
shift and spatial bins was motivated by the limited num- 
ber of independent realisations of our fiducial cosmology, 
which limits the dimensionality of the data vector. 

We predicted the cosmic shear signal on a grid 
in our paramete r spac e using the transfer function of 
lEisenstein fc HvJ (1 19981) and the non-linear power spec- 
trum of Peacock fc Doddsl |l996f) . Figure [5] shows a com- 
parison of constraints obtained using the cosmic shear 
2PCF and the he M and S peak statistics. The confi- 
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Figure 5. Comparison of the cosmic-shear tomography con- 
fidence intervals (orange/red) with the full peak information 
making use of the combined M and S statistics (blue). Shown 
in green with black outlines is the combination of cosmic-shear 
and peak statistics. Again, 1- and 2cr intervals are shown in 
all cases. 



dence region of the full peak statistics £ p is well aligned 
with the confidence region of cosmic shear tomogra- 
phy. However, as Table [2] shows, the combined statis- 
tics £ = (^p'^cs ) still gives significantly improved con- 
straints; it has a FoM that is about 40% larger than that 
of cosmic shear tomography alone. 



3.4 Stability of the constraints 



As mentioned in Sect. 12.41 the estimated inverse covari- 
ance becomes singular if d > n — 1, and the length of 
the data vectors is consequentially limited by the number 
of ray-tracing simulations at the fiducial cosmology. The 
combined statistics vector £ is 75-dimensional, which is 
comparable to the number of simulations n = 175. Even 
though we can obtain an unbiased estimate of the inverse 
covariance, this estimate is potentially very noisy. 

We estimate the effect of noise due to the finite num- 
ber simulations on the constraints derived in the pre- 
vious section by we creating 1000 bootstrap-like sam- 




Figure 6. Variation of the constraints when bootstrapping 
the covariance matrix. The white contour line shows the 95% 
confidence interval of the combined peak and cosmic shear 
statistics from Fig. \E\ Shaded areas mark the regions inside 
the bootstrapped 95% confidence interval for x percent of the 
resamplings, where x is 100 (black), 99 (red), cyan(95), and 
blue (68). 



pies from our set of simulations for the fiducial cosmol- 
ogy. Each sample is constructed by randomly drawing 
n = 175 ray-tracing realisations with replacement. The 
covariance matrix is estimated for each of the samples. 
In doing so we kept track how many independent simu- 
lations n entered the estimation of the inverse covariance 
matrix in Eq. (|12|l . At the same time we ensured that we 
had enough independent simulations to estimate C^ 1 . Fi- 
nally, we computed the corresponding confidence regions 
as described before. While not statistically rigorous, this 
scheme nevertheless illustrates that the confidence inter- 
vals are stable with respect to the set of simulations used. 
Figure [6] shows the variation of the area inside this confi- 
dence interval and illustrates that the dependence on the 
simulations entering the computation of the covariance 
matrix is small compared to the size of the confidence 
region. 

We have also compared the confidence regions ob- 
tained using the fitting functions of Appendix [A] for the 
function M with the constraints computed using the in- 
terpolation method described in Section 13.11 While in- 
terpolating between simulations for different cosmologies 
is noisy, the results using the fitting function might be 
affected by accuracy and choice of the fitting functions. 
However, we do not find significant differences between 
the two methods. 
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4 SUMMARY AND DISCUSSION 

We showed in this paper that the number and proper- 
ties of peaks found in a weak lensing survey can con- 
strain cosmological parameters. This allowed us to turn 
an important limitation of weak-lensing cluster searches, 
namely their low purity and completeness, into a source 
of cosmological informatio n. We note that a similar idea 
was recently presented by IWang et all (|20Q9h . The most 
important differences between their and our works are 
our purely numeric ansatz, the inclusion of tomographic 
information, and the combination with cosmic shear to- 
mography. 

In the pilot study presented here we demonstrated 
that the peak statistics is able to provide constraints 
on fi m and as, which are competitive with those ob- 
tained from cosmic shear tomography on the same data 
set. By comparing the constraints obtained from peaks 
found in maps of aperture mass and tomography cubes 
we concluded that the tomographic redshift leads to much 
tighter constraints than working with a 2-dimensional 
M ap - or convergence map alone. Even though the lens- 
ing efficiency varies only slowly with lens redshift and the 
best estimate of a single p eak's redshift has a large sc atter 
around the true redshift (jHennawi fc Spergell 120051 ) . the 
peak redshift distribution contains valuable information. 

By looking at the SNR distribution function we 
showed that the SNR of peaks, acting as a proxy for mass, 
also provides information beyond the simple counting of 
pe aks above a detect i on or mass threshold as it was don e 
bv iTakada fc Bridle! (|2007l > and iFang fc Haimanl (|2007l ). 
We showed that peaks with a low SNR, most of which 
cannot be associated with a single massive halo, con- 
tribute significantly to the information content of the 
peak statistics. 

Finally, we investigated how the peak statistics 
compares with the standard cosmic-shear tomography 
method and whether a combination of both statistics is 
useful. We found that the peak statistics gives constraints 
on Qm and as that are competitive with those expected 
from a cosmic-shear study on the same survey fields. Tak- 
ing at face value, our FoM suggests that the peak statis- 
tics is indeed the superior method. The comparison, how- 
ever, is not entirely fair because we used only two redshift 
bins for the 2PCF tomography while the peak statistics 
made use of almost perfect redshift information, when we 
discretized the exact redshifts of the background galax- 
ies into bins of widths A z = 0.01. On the other hand, 
iHennawi fc Spergell (| 2005h found that even with very few 
redshift bins, the tomographic peak finder is able to lo- 
cate peaks along the z-axis reasonably well and signifi- 
cantly outperforms a simple M ap peak finder. We thus 
conclude that the peak statistics gives constraints that 
are at least comparable to those obtainable from more 
traditional weak lensing methods. 

We combined the peak statistics with the cosmic- 
shear 2PCF tomography, including the cross-covariance 
of both statistics measured at our fiducial cosmology. Al- 
though both methods have a very similar degeneracy, 
their combination improves our FoM by about 40%. We 
surmised that this is due to the inclusion of information 
about the non-Gaussianity of the matter density field in 



the peak statistics, which is not contained in the 2PCF . 
This result is not unexpected. ITakada fc Bridle! (|2007t ) 
already found that the combination of cluster counts 
and shear tomography, including the full covariance, im- 
proves cosmological constraints. The new information we 
add here is that the gain in information continues to 
be present when cluster counts are replaced with peak 
counts. The peak statistics, not discriminating between 
massive halos and chance projections, acts as a "very 
high-order" cosmic shear method. 

Eventually, the much more ambitious question one 
wants to answer is: What is the ideal way to extract (lens- 
ing) information from present and future cosmological 
surveys? We made no attempt to tackle this general prob- 
lem but showed that a specific higher-order statistics, the 
peak statistics we introduced in this paper, provides in- 
formation that can be combined with existing two-point 
statistics to improve constraints. 

A great advantage of our numeric approach to the 
peak statistics is that observational effects can be in- 
cluded by the simulator to an almost arbitrary degree 
of precision, given of course these are known and un- 
derstood. An obvious example would be the inclusion of 
the holes and gaps in survey fields caused in real data 
by bright stars, diffraction spikes, ghost images, satellite 
tracks, and so on. 

The disadvantage of this numeric method is its enor- 
mous computational cost. The generation of initial condi- 
tions, the iV-body simulation, ray-tracing, and the tomo- 
graphic peak finder need about 150 CPUh per simulation 
on a cluster of Itanium Montecito 1.6 Ghz CPUs. Most 
of this time is used for the iV-body simulation. Since the 
aperture mass smoothes the convergence field with a ker- 
nel of typically 2Mpc radius high spatial resolution is 
not required. In this case the computing time can be sig- 
nificantly reduced by replacing the TreePM code with a 
simple particle-mesh algorithm. But even then, and in 
the case of massive parallelisation, the wallclock time re- 
quired to simulate one point in parameter space would 
prohibit running a standard Markov Chain Monte Carlo 
meth od. Population Monte Carlo (PMC, IWraith et aD 
120091 ) allows to investigate sample points independently 
and thus facilitates the parallelisation beyond the limits 
of effective interprocess communications of a single sim- 
ulation. 

Another approach to reduce the computation time 
of Monte Carlo simulation s with jV-body sim ulations is 
the framework proposed bv lHabib etafl j2007h . They de- 
veloped a mechanism by which the parameter space of 
interest is optimally sampled with relatively few high- 
precision simulations. Predictions for untried positions 
are made from emulations rather than full simulations. 
An important feature of this method is that the error 
bounds of the emulations are constrained and that the 
uncertainties of the emulator output are taken into ac- 
count in the parameter estimation. Either method, or 
the combination of PMC with emulators with controlled 
error, bounds should make the application of the peak 
statistics to current and future surveys feasible with cur- 
rent high-performance computing hardware. Finally, we 
point out that the growth of computing power expected 
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from Moore's Law is faster than the growth in etendue 
of envisioned survey facilities. 
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APPENDIX A: FITTING FUNCTIONS FOR 
TOMOGRAPHIC PEAKS 

We find that each component of Af(f2 m , as) is well de- 
scribed by the function 

M z (n m ,a 8 ) =Ag(z) [l + (l + zf] nia%P(z,n m ,a 8 ) 

+ P3(fim,0- 8 ) 

(Al) 
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(A2) 



where g(z) is the distance rat io Dds/Ds averaged ov er the 
source redshift distribution l)Schneider et al.lll998h . The 
functions P and P3 are polynomials, where 

Pi = Pio + PnAfi m + p 12 Aa a 
Pi = P20 + P21 AO m + P22 Acr 8 
P3 = P30 + P31 Af2 m + p 32 A(T 8 
P = 1 + z(pi + 2;p 2 ) , 

with Af2 m = f2 m — f2 mo , Aug = <jg — (jg . The constant 
A, the polynomial coefficients p m n, and the exponents 
a and /3 are the free parameters of the fitting function, 
which are fitted simultaneously for all redshift bins and 
cosmological parameters. The values of the fit parameters 
depend on the choice of the signal-to- noise-cut, the ellip- 
ticity dispersion, the redshift distribution of the galaxies, 
the filter function and radius and the galaxy number den- 
sity. Quoting the best-fit values for our specific choices of 
these parameters would therefore be of very limited use. 

The SNR probability distribution p(S; Q m , erg) is well 
described by a log-normal distribution for x = \og(S 2 /10) 

1 



p(:r;f2 m ,(T8 



exp 



og(x) - n(Q m , a s )] 2 

2<7 2 (fi m , ( J 8 ) 



(A3) 



where the cosmological dependence of /j, and a are de- 
scribed by (double) power laws 



|Ll(f2m,Cr8) —AiQ^ag 1 +QmO'8 1 



(A4) 



The cosmology dependent number count of peaks is given 
by 

AT(n m , a s ) = A 3 Q^a b 8 3 + tl^a* 3 . (A5) 

Note that we do not fit the probability density function 
(|A3|) to the data, but the cumulative distribution function 
which we compute numerically. 
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